Mon. Not. R. Astron. Soc. 000, (2009) Printed 27 June 2009 (MN WT$L style file v2.2) 



The Dynamics of Neptune Trojan: I. the Inclined Orbits 



Li-Yong Zhou 1 *, Rudolf Dvorak 2 , Yi-Sui Sun 1 

1 Department of Astronomy, Nanjing University, Nanjing 210093, China 

2 Institute for Astronomy , University of Vienna, Tiirkenschanzstr. 17, A-1180 Wien, Austria 



Accepted . Received 



ABSTRACT 

The stability of Trojan type orbits around Neptune is studied. As the first part of 
our investigation, we present in this paper a global view of the stability of Trojans on 
inclined orbits. Using the frequency analysis method based on the FFT technique, we 
construct high resolution dynamical maps on the plane of initial semimajor axis ao 
versus inclination Zq. These maps show three most stable regions, with i in the range 
of (0°, 12°), (22°, 36°) and (51°, 59°) respectively, where the Trojans are most probably 
expected to be found. The similarity between the maps for the leading and trailing 
triangular Lagrange points L4 and L5 confirms the dynamical symmetry between 
these two points. By computing the power spectrum and the proper frequencies of 
the Trojan motion, we figure out the mechanisms that trigger chaos in the motion. 
The Kozai resonance found at high inclination varies the eccentricity and inclination of 
orbits, while the vs secular resonance around iq ~ 44° pumps up the eccentricity. Both 
mechanisms lead to eccentric orbits and encounters with Uranus that introduce strong 
perturbation and drive the objects away from the Trojan like orbits. This explains the 
clearance of Trojan at high inclination (> 60°) and an unstable gap around 44° on 
the dynamical map. An empirical theory is derived from the numerical results, with 
which the main secular resonances are located on the initial plane of (ao, io)- The fine 
structures in the dynamical maps can be explained by these secular resonances. 

Key words: Planets and satellites: individual: Neptune - Minor planets, asteroids - 
Celestial mechanics Method: miscellaneous 



1 INTRODUCTION 

1 In the restricted three-body model consisting of the Sun, a 
planet and an asteroid, the equilateral triangular Lagrange 
equilibrium points (L4 and L5) are stable for all planets 
in our Solar system. Asteroids in the vicinities of L4 and 
I/5 of a parent planet are called Trojans after the group of 
asteroids found around Jupiter's Lagrange points. Objects 
on Trojan like orbits around Mars and (temporarily) around 
Earth have been observed while the Trojan-type orbits of 
Saturn and Uranus have been proven unstable due to the 
perturbations from other planets. 

As for Neptune, the possibility of stable orbits around 
the Lagra nge points have been verified in several papers , 
e.g. ilHolman fc Wisdoml 1 19931 : IWeissman fc Levisonl 1 19971 : 
iNesvornv fc Donesl 20021 ). befor e the discovery of the first 
Neptune Trojan, 2001 QR322 (|Pittichova et al.l 120031 ). Up 
to now, 6 Neptune Trojans have been discovered!]] around 
Neptune's L4 point. We list their orbital properties in 
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1 IAU: Minor Planet Center, http://www.cfa.harvard.edu/iau/ 
lists/NeptuneTroj ans .html 



Table 1. It is suspected that there could be much more 
Trojan- type asteroids sharing the orbit with Neptune than 
with Jupiter, both in the sense of number and total mass 
jSheppard fc Truiillol 120061 ). After these discoveries, more 
papers were devoted to the issue of Neptune Trojans, fo- 
cusing either on the orbital stabi l ity and origin of spe - 
cific objects (|Brasser et all l2004bl : iLi. Zhou fc Sunl I2007D 
or on creating a global view of stable regions around the 
Lagrange points, e.g. (|Marzari. Tricarico fc Scholj l2003al : 
IDvorak et aLll2007l : iDvorak. Lhotka fc Schwar j|2008l ). Both 
the observing and the theoretical studies of Neptune Tro- 
jans could give important clues on how these objects were 
trapped into their current orbits, where and when the planet 
Neptune formed and how the orbits of the outer planets 
evolved in the early stage of the formation of the Solar sys- 
tem. Therefore, it is worth to investigate the orbital stability 
of fictitious Trojans in the whole parameter space. 

All the Trojans in Table 1 are on near-circular orbits 
(with quite small value of eccentricity) and two of them have 
high inclination values. The stability and origin of inclined 
orbit is an interesting topic (|Li. Zhou fc Sunll2007l ). As the 
first part of our investigation of the whole phase space, we 
study in this paper the orbital stability of Trojans on in- 
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Table 1. Orbits of 6 observed Neptune Trojans, given at epoch 
JD=2454800.5 with respect to the mean ecliptic and equinox at 
J2000.0. The perihelion argument ui, ascending node Q and incli- 
nation i are in degrees. 



Designation 


M 


id 


a 


i 


e 


o(AU) 


2001 QR322 


57.88 


160.8 


151.6 


1.3 


0.031 


30.302 


2004 UP10 


341.28 


358.5 


34.8 


1.4 


0.028 


30.212 


2005 TN53 


287.04 


85.7 


9.3 


25.0 


0.065 


30.179 


2005 T074 


268.10 


302.6 


169.4 


5.3 


0.052 


30.190 


2006 RJ103 


238.64 


27.1 


120.8 


8.2 


0.028 


30.077 


2007 VL305 


352.88 


215.2 


188.6 


28.1 


0.064 


30.045 



clined orbits and try to find out the possible regions where 
the potential primordial Trojans could survive up to present. 
Using the frequency analysis method, we construct dynam- 
ical maps to locate the most stable regions, and figure out 
the mechanisms that bring instability to the motion. 

Since the L$ point of Neptune is nowadays in the di- 
rection of the Galaxy center thus not suitable for asteroid 
observing, it is not astonishing to see all asteroids listed in 
Table 1 are around the L4 point. Observations show that 
there are more objects around Jupiter's L4 than the L5 
point, and such an as ymmetry between L4 and L5 was also 
reported for Neptune |Holman fe Wisdomlll993T) . The origin 
of this asymmetry is discussed too in this paper. 

This paper is organized as follows. In Section 2, we 
introduce the dynamical model and the spectral analysis 
method applied in our study. We also show that the ap- 
parent asymmetry between L4 and L5 is only an artificial 
effect from the asymmetrical selection of initial conditions. 
Section 3 presents the dynamical maps around the L4 and 
L5 points. We describe the structures seen in the maps, and 
show that the Kozai resonance and the vg secular resonance 
are primarily responsible for the interesting features of the 
maps. In Section 4, the dynamical spectra of motion are con- 
structed and a semi-analytical theory of secular resonance 
is derived. Finally, the conclusions are given in Section 5. 



2 MODEL AND METHOD 
2.1 Dynamical model 

We numerically simulated the orbital evolution of thousands 
of test particles around the Lagrange points of Neptune and 
investigated their orbital stability with a method based on 
the spectral analysis on the outcome of numerical integra- 
tions. Our dynamical model includes the Sun (with the in- 
ner planets masses added onto), four outer planets (Jupiter, 
Saturn, Uranus and Neptune) and the Trojans. The Sun 
and planets gravitationally interact among themselves and 
on the Trojans, but each Trojan is assumed massless and 
therefore has no effect on other bodies. 

Since the Lagrange points are defined in the restricted 
three-body problem, we set the initial orbital elements of 
the fictitious Trojans referring to Neptune's orbit. A Trojan- 
like asteroid shares the same orbit with the planet and they 
are in fact in a 1:1 mean motion resonance. For Neptune's 
Trojan, the critical argument of this resonance is 



where \ = lu + £1 + M is the mean longitude and the sub- 
script '8' denotes Neptune (hereafter the orbital elements of 
planets are subscripted with '5, 6, 7' and '8' for Jupiter, Sat- 
urn, Uranus and Neptune as usual). We initially set 00 = 
where a c is the center of the tadpole orbits. This center is 
near ±60° for near-circular and coplanar orbit and varies 
with the semimajor axis, eccentricity and inclinat ion of the 
Trojan. Since a c changes slightly with a and i (jNamounil 
ll999l : lNesvornv et alJl2002h . we may still set a c = ±60° ap- 
proximately. 

Referring to the previ ous analysis dNesvornv fc Donesl 
120021 : iNesvornv et alJl2003 ) and to ensure a coverage of the 
representative region around the Lagrange points, the choice 
of the initial angular variables for the fictitious Trojans is 
made as below. The argument of perihelion u)o is set as cjq — 
uig — 60° for L4 (—60° for L5), while the ascending node tto 
and the mean longitude Mo are the same as that of Neptune 
fio = Cls,Mo = Ms. After some test simulations, the initial 
semimajor axes ao of the test particles are determined to be 
evenly sampled in the range from 29.72 AU to 30.32 AU for 
L 4 and from 29.90AU to 30.50AU for L B . Note these two 
ranges of a are different from each other (see a discussion in 
the following section). 

The observed Neptune Trojans all have small eccen- 
tricities but their inclinations are widely spread, particu- 
larly, two of them have high inclinations 25° and 28°. In 
this paper, we focus on the dynamics of inclined orbits, and 
therefore we fix the initial eccentricity of test particles at 
eo = e$ = 0.00625 w (nearly circular) and vary the initial 
inclination io from 0° to 70° with an increment of 1.25°. 

In the above arrangement, the initial value of the res- 
onant argument a is always at the libration center. For a 
Trojan, the argument a generally librates around the cen- 
ter o c with a nonzero amplitude D, and simultaneously the 
semimajor axis a oscillates with an amplitude d around a 
mean value g iven by the semimajor a xis of Neptune. It has 
been shown !|Erdilll988l ; iMilanil Il993h that the values of d 
and D are constant in the simplest approximation and they 
are related to each other by: 

d 



3/xag + 



0.373, 



(2) 



a = A — A« 



(1) 



where [i is the mass of Neptune measured in the solar mass, d 
is measured in AU and D in radians. For this reason, having 
different ao values, we do not need to test different initial 
value of a. 

Given the proper initial condit ions, the systems are 
then integrated with a Lie-integrator (jrlanslmeier fc Dvorak! 
Il984h . In some cases, we use also the hybr id symplectic in - 
tegrator in the Mercury6 software package (|Chamberdll~999h 
to compare the results and to simulate the orbital evolution 
for the lifespan of the solar system (4.5 Gyr). 



2.2 Spectral analysis 

During the integration of the system, the orbits of the fic- 
titious Trojans are recorded for further analyzing. The Lie- 
integrator outputs the orbits at every equal time interval so 
that we can apply a Fast Fourier Transform (FFT) method 
to the orbital variables. The power spectra of these variables 
then can be used to give information about the dynamical 
behavior of the orbits. On one hand, the system must be 
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integrated to a time long enough to reveal the secular be- 
havior; on the other hand, if we record the orbits in such a 
long time span, the data to be analyzed would be by far too 
large. 

The amount of the output data can be reduced by in- 
creasing the time interval at which the output is recorded. 
But the risk of losing valuable information and/or intro- 
ducing artificial features arises from the large sampling pe- 
riod. As a compromise, a low-pass digital filter is intro- 
duced to remove the short-period terms, which generally 
are not crucial and can be cut off by the average pro- 
cedure ijCarpino. Milani fc Nobilil Il987l ). Then we collect 
the records from the filtered output at every tens-of-times 
sampling intervals thus the amount of data is largely re- 
duced. A digital filter s upplied by T. Michtchenko is ap- 
plied in this paper. See (|Michtchenko fc Ferraz-Mellol 1 19951 ; 
iMichtchenko et al.| [2002) . where the details about this filter 
and its application can be found. 

After some test runs, we carefully choose an output in- 
terval of 4yr. The output data are then smoothed through 
the filtering process and among the filtered data a sampling 
interval of A = 512 yr is adopted. This interval is shorter 
than all the secular periods in the outer solar system, such 
that all information on the long-term effects is preserved. 
Finally, the systems are integrated to ~ 3.4 x 10 7 yr and 
totally N = 65, 536 (= 2 16 ) lines of the filtered orbital vari- 
ables including a,e,i,a etc are recorded for further ana- 
lyzing. The sampling interval determines the Nyquist fre- 
quency /Nyq = 2k = 9-^66 x 10 _4 yr _1 , while the fre- 
quency resolution calculated from this data through an FFT 
is /ros = = 2.980 x 10~ 8 yr -1 . We have tried smaller out- 
put intervals that can cover short period effects such as the 
quasi 5:2 mean motion resonance between Jupiter and Sat- 
urn (the Great Inequality, with a period of ~ 880 yr) but did 
not find any crucial effects with periods shorter than 10 3 yr 
in the motion of Neptune Trojans. 

An FFT to the output data provides us valuable in- 
formation about the orbital dynamics, and it allows us to 
determine the regularity of an orbit. Generally, a regular 
orbit is characterized by quasi-periodic terms, so that the 
power spectra of the orbital elements of this orbit are dom- 
inated by a countable (small) number of frequency com- 
ponents. On the contrary, a chaotic orbit is not quasi- 
periodic and the independent frequencies of the motion 
change with time, such that the power spectrum consists 
of broadband components and is characterized by strong 
noise. To tell the difference between the spectra of regu- 
lar and chaotic motion, we can count the spectral peaks 
above a given noise level, and the number obtained is de- 
fined as the spectral number (SN hereafter). Obviously, a 
small SN indicates a regular motion while a large one 
implies a chaotic motion. The SN has been successfully 
appli ed to indicate the regularity of the main-belt aster- 
oids ijMichtchenko fc Ferraz-M ello 1995; Michtchenk o et al] 
2002l')and extra solar planets dFerraz-Mello et al.l 2005; 



Michtchenko. Beauee fc Ferraz-Melldl2008l ). We use this in- 
dicator to construct the dynamical map in Section 3. 

2.3 Z/4 versus L5 

Since long ago, the asymmetrical distribution of stable re- 
gions around Neptune's L4 and L5 points puzzles. For exam- 
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Figure 1. Upper panel: The temporal variation of the osculating 
semimajor axes of Neptune and of the Trojan at the resonance 
center. Lower panel: the power spectrum of the Trojan's semima- 
jor axis in the above panel. 



pie, a shift in the semimajor axis of stable regions around 
L4, and L5 was found when test partic les in the outer so- 
lar sy stem were integrated t o 20Myr dHolman fc Wisdoml 
1993). Nesvorny and Dones iNesvornv fc Pones! (|2002h ar- 
gued that the asymmetrical distributions of stable region 
was nothing more than an artificial effects rising from selec- 
tions of initial conditions. When we set the initial conditions 
for test particles referring to the osculating orbit of Nep- 
tune, generally they are symmetrically distributed around 
the Z/4 and L5 points. It is worth emphasizing that this 
symmetry is only valid in the frame consisting of the Sun, 
Neptune and the asteroid. But the real system evolves sym- 
metrically with respect to the bary center, which is not ex- 
actly in the Sun. In this sense, the initial conditions around 
Li and L5 are not symmetric to each other any longer. 
Such an argument was also prop osed by Marzari et al. 
M arzari. Tricarico fc Scholll ll2003ah. Mo re recently, Dvorak 
et al. iDvorak. Lhotka fc Schwarzl ( 20081 ) showed by sophis- 
ticated numerical experiments that a proper selection of ini- 
tial conditions can remove the difference between L4 and L5. 
Nevertheless, we would like to show here another evidence: 
we construct separately two dynamical maps for Trojans 
around the L4 and L5 (see Figs. 2 and 3). Judging from the 
appearances, they are symmetrical to each other, only ex- 
cept a shift in the semimajor axis, which we analyze below. 

Given the initial orbital elements of 101 test Trojans as 
eo = es,io = «8,<^o = + 60°,Slo = Q.s,Mq = Ms and 
ao varying from 29.9 AU to 30.5 AU, we integrate the sys- 



tem up to 3.4 x 10 4 yr (the typical libration period o f a is 
~ 9x 10 3 yr, see Chap. 3.8 of (|Murrav fc Dermottlll999h ) and 
calculate the amplitude Act of the librating resonant argu- 
ment a for each Trojan. The value of ao at which Act meets 
the minimum, assigned a™ 1 ", is regarded as the center of 
the tadpole orbits around L4 point. At a specific moment t, 
the osculating orbits of planets lead to one value of a nln (i). 
At next moment t' , the planets have evolved to new oscu- 
lating orbits. We restart the above-mentioned calculations 
with planets on the new osculating orbits, but the fictitious 
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Figure 2. The dynamical map around L4 point. The color indi- 
cates the spectral number. Green indicates regular motion, while 
red is on the edge of chaotic motion. From our experiences (see 
Section 3.2), those orbits with SN > 50 are roughly not to be ex- 
pected to survive in the lifespan of the solar system. The libration 
amplitude Act of the resonant angle is also mapped by contour 
curves. Among them, the red contour curve represents Act = 10° 
while the thick black one is for Act = 60° . 

Trojans still on the same initial orbits as before. A new value 
of the central semimajor axis a™ m (t') is obtained. In such 
a way, we finally get two series of semimajor axes, one is 
for Neptune's osculating orbit and the other for the corre- 
sponding Trojan at the resonance center. We show the time 
variations of them in Fig. 1. 

Applying an FFT to the time series of the central a™ m 
in the upper panel, we get the power spectrum in the lower 
panel of Fig. 1. The frequencies of the peaks tell us the 
mechanism causing the variation of a™ ln . The highest four 
peaks are at f 1 = 0.0781, f = 0.0273, f = 0.162 and 
f = 0.00549 (27r/yr). Denoting the mean motion (orbital 
frequency) of planets by fs,fe, JV and /s , a simple calcu- 
lation reveals that f 1 = / 5 - fs,f 2 = /s - /s,/ 3 = 2/ 5 
and / 4 = fy — fs, that is, these frequencies are either the 
synodic frequencies (J 1 ' 2 ' 4 ) or the harmonics of orbital 
frequencies in the outer solar system. This fact demonstrates 
again that the asymmetrical "shift" in the semimajor axis 
of the center of the tadpole orbits is only due to the initial 
orbital configurations of planets. Because of its largest mass, 
Jupiter plays the most distinct role in causing this variation. 



3 RESULTS 

3.1 Dynamical map 

We present our investigations on the dynamics of the in- 
clined Trojans in this part. As mentioned above, the value 
of the libration center <r c changes only slightly with inclina- 
tion, and therefore we may fix its value at 60° for L4 and 
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Figure 3. The same as Fig. 2 but for L5 point. 



—60° for L5 when we study the dependence of stability on 
the orbital inclination. 

Adopting the spectral number as an indicator of the 
regularity of orbits, we construct dynamical maps using 5757 
orbits starting from a 101 x 57 grid on the (ao>*o) plane 
and integrated for 34Myr. The power spectrum of cos cr for 
each orbit is calculated and the number of peaks, which 
are over one percent of the highest peak, is defined as the 
SN. To limit the number into a manageable range, an SN is 
forced to be 100 if it was originally larger than that. We also 
exclude orbits that escape from the 1:1 resonance region. 
A simple criterion is applied: if the averaged value of the 
semimajor axis a of a test particle does not satisfy 29.9 AU < 
a < 30.5 AU, it is regarded as escaped from the resonance, 
and an SN of 110 is assigned to the orbit. We show the 
dynamical maps in Fig. 2 for L4 and in Fig. 3 for L5 . Because 
all orbits with inclination higher than 61° can not survive 
in the Trojan-like orbit, we show in Figs. 2 and 3 only for 
orbits with i G [0°, 63.75°]. 

The colour in the dynamical map indicates the SN. In 
the green region where SN is relatively small, the motion is 
dominated by a few dominating frequencies and thus it is 
more regular; but in the blue and red region, the spectrum 
of cos a is characterized by strong noise and thus the motion 
is chaotic; the white color with SN = 110 indicates escaped 
orbits. 

The dynamical maps in Figs. 2 and 3 do not show any 
asymmetry between L4 and L5. The only distinguishable 
difference between them is a shift in the semimajor axis, 
which is due to the setting of initial conditions, as described 
in Section 2. Hereafter in this paper, we will discuss only the 
Trojans around the L5 point (Fig. 3). We believe that the 
results can be mirrored to the L4 point. 

Judging from the dynamical maps, the most regular or- 
bits mainly exist in three regions in inclination: A: 0° — 12°, 
B: 22° — 36°, and C: 51° — 59°, approximately. The region 
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A and B are connected by a less regular region, but re- 
gion C is apparently separated from others by an unsta- 
ble gap at io ~ 44°. Note the inclination values of the ob- 
served Trojans in Table 1 in fact reside in region A and 

B. Except for the white-colored region, our results do not 
repel the possibility of Trojans outside regions A, B and 

C, but we argue that the probability of finding Trojans 
there is lower. Particularly, region C at high inclination 
is isolated from regions A and B, therefore any Trojans 
found in region C should have been always there, i.e. they 
must be primordial. The upper limit of region C is con- 
sistent with the value (61.5°) given in a rest ricted three- 
body model l|Brasser. Heggie fc Mi kkola 2004a!'). The upper 
limit of region B is very close to the value (35°) given in 
(|Nesvornv fc Donesl 120021 ). where the authors argued that 
no Trojan with higher inclination could survive. By plot- 
ting "diffusion maps" f or initial inclination ip at 0° , 10° , 20° 
and 30°, Marzari et al. iMarzari. Tricarico fc Scholll (1200381 ) 
found two "high stability regions" locating respectively at 
io — 0° and io = 30°. They have used the slice of phase 
space at definite inclinations, and the two stable slices re- 
ported by them are well inside region A and region B. 

The coplanar orbits are not necessarily more regular 
than the inclined ones. In fact, there are two less regu- 
lar "holes" at low inclination io < 5° and around ao ~ 
30.13, 30.31 AU in Fig. 3. The central area, in which the li- 
bration amplitude is small, is likewise not necessarily more 
safe for Trojans. Different secular resonances are crowded in 
this area and their overlappings may introduce chaos. 

We also calculated the libration amplitude of the reso- 
nant argument a. Here by "amplitude" we mean the differ- 
ence between the maximum and minimum values of a during 
our integration time, Act = CT max — CT m i n . The contour of Act 
is also plott ed in Figs. 2 and 3. It has been shown in previous 
litera tures l|Weissman fc Levisonl 1 19971 ; INesvornv fc Donesl 
l2002h that Act should not exceed 60° - 70° for stable or- 
bits in the resonance. In fact the green (regular) region in 
Figs. 2,3 is shaped by the contour curve of Act = 60°. But for 
small inclination io < 10° there are some regular orbits hav- 
ing Act ~ 70°. The most regular orbits in region A, B and 
C have libration amplitude of 20° — 60° . Submit Act = 60° 
(D — it/3) into Eq.(2), we get the oscillating amplitudes of 
semimajor axis d ~ 0.39 AU. It is a little larger than the 
value obtained from Fig. 3 (~ 0.31 AU), which is shrunk by 
perturbations from other planets besides Neptune. The li- 
bration amplitude up to 60° and even to 70° indicates also 
that all the Trojans which survive in our integrations are on 
the tadpole orbits. Trojans initially on horseshoe orbits in 
our simulations typically escape from the resonance before 
10 Myr. 

The maximal eccentricity of an orbit experienced 
during the evolution, e max , has bee n used as an indi- 
cator of the chaoticity of mo tion (|Dvorak et al.l 120071 ; 
iDvorak. Lhotka fc S chwarz 2008]). We see from such an e max 
map in Fig. 4 that most of the orbits in the 1:1 mean mo- 
tion resonance maintain their small eccentricity during the 
evolution. In those areas corresponding to green regions in 
Fig. 3, the eccentricity is smaller than 0.05. An interesting 
feature visible in Fig. 4 is a strip of relatively larger eccen- 
tricity crossing through the map at io — 11° — 16°. In this 
strip, e ma x ~ 0.05 while in its vicinity we have e max ~ 0.03. 

There are still some other fine structures buried in the 
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Figure 4. The maximal eccentricity e max of orbits experienced 
in the evolution. 

dynamical map in Fig. 3 (and also, similar in Fig. 2). For 
example, an "arc" of less regular motion than its neighbor- 
hood, with ends at (a = 30.04AU,i = 0°), (30.38AU,0°) 
and crossing (30.20 AU, 22°), can be easily recognized in 
Fig. 3. This notable feature is also characterized by a con- 
striction of the Act contour curves, i.e. orbits initialized on 
this arc have larger libration amplitudes. A closer look at 
this region reveals that this unstable arc is encompassed ex- 
teriorly by a stable arc and an unstable but shorter arc. 
Another "arc structure" can be seen extending from (ao = 
30.02 AU,i = 20°) to (30.08 AU, 30°), which is also accom- 
panied by a slight deformation of the Act contours. Last 
but not least, both the dynamical maps in Figs. 2 and 3 are 
symmetric with respect to the central line (ao = 30.02 AU 
in Fig. 2 and ao = 30.20 AU for Fig. 3), where locates the 
center of the tadpole orbit. This symmetry is apparent 
but not exact. The tiny deviation from the symmetry is 
well-known and has been described, e.g. in Chap. 3.9 of 
l|Murrav fc Dermott|[l999l '). 

So far we have depicted the orbital stabilities of Tro- 
jans on inclined orbits with their initial eccentricities all re- 
stricted to a small value (nearly circular orbits). Certainly, 
the stability of a Trojan's orbit depends on its eccentric- 
ity. In fact, the stable region in eccentric i ty is very lim- 
ited. Nesvorny & Pones INesvornv fc Donesl (|2002h reported 
that for inclination i = 0° all the "low-LCE" (stable) tad- 
pole orbits have e < 0.1, and they also found an orbit of 
e = 0.07, i — 25° surviving over 4 Gyr in their numeri- 
cal inte grations. Marzari et al. IMarzari. Tricarico fc Scholll 
l|2003ah showed that in the "high stability regions" on the 
slices of io = 0° and io = 30°, the proper eccentricities of 
Trojan orbits are small, e p < 0.1 and e p < 0.15 respectively. 

We also calculated the dynamical maps on the (ao,eo) 
plane with specific inclination values. Our results confirm 
the conclusion that the eccentricities of stable orbits are re- 
stricted to small values. As examples, we show in Fig. 5 two 
cases with io = 5° and io = 55°, correspondingly inside the 
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Figure 5. The dynamical maps on the (a, e)-plane for specific 
initial inclinations. The upper panel is for io = 5° and the lower 
for io = 55°. The color codes are the same as in Fig. 2 and 3. The 
contour curves represent the libration amplitude of the resonant 
argument ct. From the inside out are contours for Act = 10° to 
Act = 70° with an increment of 10° . 



stable region A and C. For io = 5°, an orbit with initial 
eccentricity as large as eo ~ 0.085 may still be among the 
most stable orbits and its libration amplitude may reach 
Act ~ 70°. For io = 55°, however, the most stable region 
extends only below eo ~ 0.02 and the libration amplitude 
Act < 60°. Our preliminary calculations for other slices at 
different inclination values show that some test particles 
with eo ~ 0.25 at io = 35° may retain on the Trojan orbits 
in our 34Myr integrations, but the SN indicator indicates 
that the eccentricity of the most stable orbit in region A, B 
and C, approximately, should not exceed 0.10, 0.12 and 0.03 
respectively. 

Only two slices have been shown here and apparently 
in Fig. 5 there are fine structures that bear plenty informa- 
tion about the orbital dynamics. Toward a global view of 
the stability of eccentric orbits, particularly, of the mecha- 
nisms behind, we need more slices at different inclinations 
and a through analysis on the orbital dynamics. This work is 
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Figure 6. The lifetime (in logarithm) of orbits from the 4.5 Gyr 
integrations. These orbits are initialized either on the horizontal 
line (upper panel) or the vertical line (lower panel) on the (ao,io) 
plane (see text for an explanation). Their spectral numbers are 
calculated from the short-term (3.4 X 10 7 yr) integration, and in- 
dicated by color. Note those orbits with SN = 110 (see text) is 
colored grey in this figure. 



undergoing and we would like to leave this to a separated pa- 
per. From below in this paper, we will focus on the dynamics 
of the inclined orbits, i.e., we will analyze the mechanisms 
causing the structures in the dynamical maps of Figs. 2 and 
3. 



3.2 Spectral number and long-term stability 

Before continuing to analyze the mechanisms that por- 
trait the features of the dynamical maps, we perform some 
long-ter m orbital integrat ions using the Mercury6 integrator 
package ((Chambers 1999) to compare them with the results 
from our Lie-integrator and to check the reliability of the 
indicator i.e. the spectral number. 

We arbitrarily select initial Trojan orbits from two lines 
on the (ao,io) plane. One is a horizontal line with io = 10° 
and the other one is the vertical line with ao = 30.098 AU. 
Hundreds of orbits initialized on the lines are integrated up 
to the solar system age, 4.5 Gyr. An object is discarded if 
its semimajor axis is larger than 60 AU. A Trojan can attain 
such a wide orbit through close encounters with the Sun 
and/or planets. We do not check the time when a Trojan 
leaves the 1:1 resonance, but generally, an object will be 
scattered far away soon after it loses the protection of the 
1:1 resonance. 

The results are shown in Fig. 6. We note the correlation 
between the lifetime and the SN. Those orbits surviving the 
whole integration are also those with small SN while orbits 
with relatively larger SN escape from the resonance before 
the integration ends. Particularly, in the lower panel we see 
four orbits with i = 13°, 14°, 15°, 16° lose their stability 
after 4.2, 0.24, 1.1 and 2.4 Gyr respectively. These points 
are on the blue arc in Fig. 3, i.e., the chaotic property of 
these orbits has been correctly predicted by the relatively 
large SN calculated from our short-term integration. 

We should point out that several orbits with small SN 
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Figure 7. An example of unstable orbits (ao = 30.218 AU, io = 
61.25°) trapped in the Kozai resonance. From top down, 5 panels 
are for semimajor axis a (in AU), eccentricity e, inclination i, 
argument of perihelion u and the resonant argument a. Angular 
variables are in degrees. 

(green ones in Fig. 6) can not sustain the solar system age. 
They are all at the border of stable region, and escape only 
after very long evolution (~ 10 9 yr) . Perhaps the instability 
is introduced through very slow chaotic diffusion, which the 
SN fails to detect. Nevertheless, one can see that globally 
the SN is still a successful indicator of orbital stability. 

Another conclusion can be made from Fig. 6, namely 
Trojans with high inclinations io > 50° can survive the solar 
system age. io = 35° is not the upper limit in inclination for 
potentially stable Trojans, as one also can see from Fig. 2 
and Fig. 3. 

3.3 Kozai mechanism and vs resonance 

The most distinguishable features in the dynamical map 
are two white (unstable) regions: the high inclination re- 
gion with io > 61° and an unstable gap at io ~ 44°. Orbits 
in these two regions are strongly chaotic and they cannot 
survive in the resonance even in our short-term integration. 
There must be some strong mechanisms driving them out. 

For orb i ts with high inclination, the Kozai resonance 
|Kozailll962l ; iKinoshita fc Nakaill2007r i is acting as the re- 
sponsible mechanism. In a Kozai resonance the perihelion ar- 
gument w librates while the eccentricity and inclination un- 
dergo variations such that the quantity Hk = y/i ~ e 2 cosi 
remains constant. When the inclination of a Trojan de- 
creases its eccentricity increases, as a result, it will cross 
Uranus' orbit and the probability of close encounter with 
Uranus is enhanced significantly. 

We checked the evolution of some orbits with io > 61° 
and found that the instability is really due to the Kozai 
resonance. A typical orbit is shown in Fig. 7 We see clearly 
that u> enters a small-amplitude libration state at 0.5 Myr 
after a transient period. From then on u is nearly constant, 
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Figure 8. The evolution of a typical orbit (ao = 30.218 AU, 
io = 43.75°) initialized in the unstable gap. From top down, 6 
panels show the temporal evolution of semimajor axis a (in AU), 
eccentricity e, inclination i, argument of perihelion ui, 1:1 resonant 
argument a and the difference between the perihelion longitudes 
of the Trojan and of Neptune. All angular variables are given in 
degrees. 



e increases as i decreases till 4.6 Myr when e reaches a value 
of 0.358. The perihelion distance is now q = a(l — e) w 
30.20 x (1 - 0.358) = 19.39 AU, which is in the range of 
Uranus distance to the Sun. A close encounter with Uranus 
then destroys the 1:1 mean motion resonance between the 
Trojan and Neptune, as shown by the deviation of a from 
libration in the bottom panel of Fig. 7. After that the object 
is finally scattered out from the solar system through a series 
of encounters with planets. 

Generally, the Kozai mechanism protects an object from 
close encounters with the perturber (usually a planet) ei- 
ther for orbits at high inclination where u> librates around 
90° or 270° l|Thomas fc Morbidellil Il996l ). or for orbits 
at low inclination wher e u> oscillates around 0° or 180° 
(|Michel fc Thom as 1996). But in our case of Neptune Tro- 
jans, as seen in the third panel of Fig. 7, ui librates around 
322° with an amplitude smaller than 5°. To our best knowl- 
edge, there is no such an "asymmetrical" libration center of 
w being reported before. However, in the case of Neptune 
Trojan, there are at least two facts that must be taken into 
account in the analysis of the Kozai resonance. First, the 
object is in a 1:1 mean motion resonance with Neptune, and 
second, there is more than one perturber on the asteroid's 
motion. Thus we believe the Kozai resonance here is more 
complicated and it deserves a careful investigation in future. 



8 Zhou, Dvorak & Sun 



As for the unstable gap at io ~ 44° in the dynamical 
map, a close look at the orbits reveals that it is due to the 
apsidal secular resonance iyg ■ A vs secular resonance happens 
when the precession rate of the Trojan's perihelion equals 
the fundamental frequency gg of the solar system, which is 
mainly related to the apsidal precession of Neptune. Roughly 
speaking, in a vg resonance, the Trojan's perihelion precesses 
at almost the same rate as Neptune's perihelion does, and 
the difference between the longitudes of perihelion zu — vjg 
oscillates around a constant value with a definite amplitude. 
We illustrate in Fig. 8 the orbital evolution of a typical Tro- 
jan initialized in the unstable gap to show the effects of the 
vg resonance. 

As shown in Fig. 8, a librates with a small amplitude 
around the Lagrange point L5 with a G (—65°, —52°) before 
f .57 x f 7 yr. Thanks to the protection of the 1:1 resonance, 
the evolutions of other orbital elements during this period 
are regular, e.g. a is nearly constant and i varies around 
~ 44° with an amplitude of only ~ 4°. But there is one 
exception, the eccentricity e is increasing during this period 
and it reaches e = 0.355 at T = 1.57 x 10 7 yr. We know that 
the secular resonance related to the perihelion preces sion 
may drive the eccentricity up |Murrav fc Dermottil 19991 ) . In 
fact, the vg secular resonance, characterized by a libration of 
vj — vjg as shown in the bottom panel of Fig. 8, is responsible 
for this eccentricity increasing. 

Again, the high eccentricity, this time owing to the 
vg resonance, reduces the Trojan's perihelion distance, and 
then the object is driven out by the strong perturbations 
during close encounters with Uranus. We may note that af- 
ter leaving the 1:1 mean motion resonance and before being 
scattered far away, the object temporally experiences the 
Kozai resonance again from 16 Myr to 20 Myr with uj oscil- 
lating, around 180° this time. 

3.4 Three-body resonance 

So far we discussed two mechanisms that affect the secular 
behavior of Trojan orbits, i.e. the Kozai resonance and the 
vg resonance. Apart from the unstable regions due to these 
two mechanisms, some other chaotic subregions embedded 
in the dynamical map can be seen, e.g. the blue (chaotic) 
arc structures mentioned in Section 3.1. The possible mech- 
anisms behind these structures are not so obvious. However, 
since the orbital periods of Neptune and Uranus are very 
close to a 2:1 commensurabilit y (165 yr vs 84 yr), the three - 
body mean motion resonance l|Nesornv fc Morbidellilll998l ) 
between the Trojan and them is a good guess. 

For a Neptune Trojan, we expect that a three-body res- 
onance may happen when A + As — A7 ~ 0. The resonant 
angles associated with this resonance would be any combi- 
nation of the form: 

A + As — A7 + Izu + I7ZU7 + Igzug, (3) 

where l,h,lg are integers satisfying the d'Alembert rule: 
I + If + lg = — 1. Different resonances with different com- 
binations of I, I7 and lg may be responsible for the "multi- 
plet" arc structures in the dynamical map. By testing the 
behavior of orbits, we find that several (but not all) Trojans 
initialized in the most distinct chaotic arc may be related 
to the l,l7,lg combination of 2, -3, 0. For these orbits the 
angle A + As — A7 + 2zu — 3zu7 typically varies very slowly. 



Table 2. The fundamental secular frequencies in the 
outer solar system. The pe riod values are taken from 
jNobili, Milani fc Carp"inolll98g|) . from which the frequency val- 
ues are computed. The periods are given in years and frequen- 
cies m lCT 7 27r/yr. Since the resolution in frequency is 2.980 X 
10 -8 2n/yr in our FFT procedure, we keep two places of decimals. 





Period 


Freq. 




Period 


Freq. 


95 


304,400.48 


32.85 


S5 


129,550,000. 


0.08 


96 


45,883.37 


217.94 


S6 


49,193.46 


203.28 


97 


419,858.29 


23.82 


S7 


433,059.42 


23.09 


96 


1,926,991.9 


5.19 


SB 


1,871,442.70 


5.34 



This is an "order 5" resonance. But in fact, as we will see 
in next section, those arcs probably are not caused by the 
three-body resonance. 

It is impossible and unreasonable to check all the com- 
binations of the integers l,lr,l$ for all orbits inside an in- 
teresting area. To obtain a global understanding of motion 
on the initial plane (ao,io), we turn to a frequency analysis 
method in next section. 



4 FREQUENCY ANALYSIS 

We have integrated thousands of Trojan orbits to construct 
the dynamical map and the filtered output from the integra- 
tions contains valuable information about the global view of 
the motion on the initial plane. We show in this part our fre- 
quency analysis on three basic variables in the motion: the 
resonant argument in the form of cos a and the non-singular 
variables k and q, which are related to the eccentricity e and 
inclination i through the relations: 

k = ecosru; g = sinicosf2. (4) 

The spectra of these variables give information about the 
most important rates of the resonant argument, the peri- 
helion longitude and the nodal longitude. We denote these 
three proper frequencies by fa,g and s hereafter. 

4.1 Dynamical spectrum 

Typically, a power spectrum of cos a, k or q is very infor- 
mative but simultaneously very complicated as well. It is a 
complicated composition of peaks at all the forced frequen- 
cies, free frequencies, their harmonics and their combina- 
tions. Therefore a direct look at a specific power spectrum 
may not be very helpful. However, when we exam the contin- 
uous change of the spectrum with a parameter, some impor- 
tant features emerge. This continuous change of a spectrum, 
defined as the dynamical spectrum, is calculated in this pa- 
per. 

We first list the fundamental frequencies 
jNobili. Milani fc Carp"inolll989h in the outer solar system 
in Table 2. One frequency not listed but important in 
the dynamics of Neptune Trojan is the frequency of the 
quasi 2:1 mean motion resonance between Neptune and 
Uranus, /2N:iu- The value determined from our calculation 
is /2N:iu = 2.3606 X 10 -4 2-7r/yr, corresponding to a period 
of 4236 yr. 

For each orbit initialized on a horizontal or vertical 
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Figure 9. Dynamical spectra of the apsidal variable k. The top 
panel is for initial orbits on the vertical line ao = 30.098 AU 
on the (ao,eo) plane, and the bottom panel is for orbits on the 
horizontal line io = 15°. The frequency of the highest peak is 
denoted by open squares, and the frequencies of the 2nd, 3rd and 
4th highest peaks are open circles, open triangles and solid cir- 
cles respectively. Other frequencies are given in different type of 
symbols, but not labeled. The dashed lines give the values of the 
fundamental frequencies in the outer solar systems 35 , 56 ; 97 ; 98 
and the frequency of the quasi 2:1 mean motion resonance be- 
tween Neptune and Uranus /2N:1U- The thick solid curves are 
the numerical fit of the proper frequency of k given in Eq. (5) and 
the thin solid curve stands /2N:iu ~ 2/ CT from the numerical fit. 



straight line on the (ao,io) plane, the FFT is performed on 
the filtered output of k (or q, cos a) and a power spectrum is 
obtained. The leading terms in the spectrum are picked out 
for further analysis. Here by "leading terms" we mean those 
frequencies at which the peaks in the power spectrum are 
among the highest ones. We show first in Fig. 9 the variation 
of the leading frequencies of k along the specific lines on the 
(ao,*o) plane. The strength of the peak is not specified in 
number but denoted by its order in the sequence of peaks' 
heights. 

As shown in Fig. 9, the motion is complicated with the 
spectrum typically characterized by a composition of many 
terms at different frequencies. The peaks at the fundamen- 
tal frequencies are the forced terms, and their frequencies 
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Figure 10. Dynamical spectra of the nodal variable p. The top 
panel is for initial orbits on the vertical line ao = 30.152 AU, and 
the bottom panel is for orbits on the horizontal line io = 15°. 
The dashed lines denote the fundamental nodal frequencies sq,sy 
and sg. The frequency related to the nodal precession of Jupiter, 
S5 ~ 8 X 10~ 9 27r/yr, is beyond the scope of this figure and can 
not be correctly reflected in our integration of 34 Myr. The thick 
solid curves stand the numerical fit of the proper frequency s in 
Eq.(6) and the thin solid curve in the top panel indicates the 
frequency of S7 — 2s. 

do not change with the parameters (io in the top panel and 
ao in the bottom panel) . The proper frequency can be easily 
recognized because it changes continuously with the param- 
eter. 

The proper frequency of the apsidal precession g shown 
in the top panel of Fig. 9 is always below 35, g§ and g-j, but it 
meets gg at io ~ 43.5° where the vg secular resonance (g — 
gg) is located. Around this resonance and also at the high 
inclination end, the motion is chaotic, as reflected by the 
discontinuity at the frequency evolution and the scattered 
values of frequencies. 

As described above, with the help of dynamical spec- 
trum, we can locate the position of the vg secular resonance. 
In fact, the dynamical spectrum contains more information. 
For example, a frequency at low inclination can be seen ap- 
proaching g§ with increasing io and meeting <?6 finally at 
io ~ 18° , as indicated by an arrow in the top panel. Checking 
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the dynamical map, we see this vertical line ao — 30.098 A 
crosses the "arc structure" from io = 14° to 20°. Arouri 
io — 18° is a regular island surrounded by two chaotic sej 
ments at 15° and 20°. Reasonably we suspect that the ai 
structure arises due to this frequency crossing. A careful ca 
culation tells us that the varying frequency seen in Fig. 9 
probably /sn : iu - 2/ CT . 

The proper frequency g is not necessarily the dom 
nant frequency in the motion. As we can see in Fig. 9, tl 
forced term associated with Uranus dominates the motic 
at low inclination where the highest peak (denoted by ope 
squares) in the spectrum is the forced term at frequent 
of gj, and the forced term associated with Neptune taki 
over the control at higher inclination (io > 17°). Only i 
the range 13° ^ io ^ 16°, when the value of the propi 
frequency is far from the values of both g-j and g%, the Tri 
jan's motion is dominated by its own proper frequency. Th 
explains the fact that the e max in Fig. 4 reaches its max 
mum around io ~ 14° , because otherwise the eccentricity 
suppressed by the forced oscillations. 

Similar dynamical spectra for the nodal variable q ai 
presented in Fig. 10, and a similar analysis on the dynam 
cal spectrum can be done too. The proper frequency of the 
nodal precession s is far away from any fundamental fre- 
quencies except the sg. Although it is very close to s% at 
low inclination, the effect of secular resonance i>i$ (s — sg) 
seems very weak compared to the or even some secondary 
resonances, as we will discuss in next subsection. 

The nodal precession of the Trojan is less affected by 
planets. The dominant frequency in the spectrum of q is 
the proper frequency s, and the leading frequencies include 
s, 2s, sr, S7 + 2s and S7 — 2s, as shown in Fig. 10. This im- 
plies that except in the region of vi$ resonance, the most 
strong perturbation to the nodal evolution of Trojan is from 
Uranus. 

Because of the small value (s < 6 x 10~ 7 27r/yr) and of 
the restricted resolution of FFT, the precision of s is poorer 
compared to g and / CT . In this sense, a longer integration 
time is needed to give more accurate results about the nodal 
resonance of Trojan motion. 

The dynamical spectrum of the resonant argument a 
can also be constructed in a similar way. We do not show 
them here but present empirical expressions of the proper 
frequencies g, s, f a and construct a resonance map on the 
initial plane (ao,io). 

4.2 Semi-analytical model 

Knowing the values of the proper frequencies g,s and 
on the initial plane (ao,io), we obtain empirical expressions 
through numerical fitting. Set x — ao — 30.219, y = sinip 
and adopt the quadratic formula used in references dMilanil 
1 1994 iMarzari. Tricarico fc Scholll l2003bft . the best fit of 
these proper frequencies from our calculations are: 

g[W' 6 27r/yr] = 1.586 + 5.614x 2 - 2.747y 2 

- 40.386a; 4 + 1.192y 4 - 15.811a;V (5) 
s[10~ 7 27r/yr] = 5.078 + 19.791x 2 - 8.218y 2 

- 54.957a; 4 + 4.192j/ 4 - 7.855a; V (6) 
/ CT [10~ 5 27r/yr] = 11.349 - 17.421a; 2 - 3.846y 2 
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Figure 11. The main secular resonances (commensurabilities) for 
Trojans around the L5 point on the (ao, io) plane. The resonances 
are labeled along the curves. Three thin curves labeled 1, 2 and 
3 represent the locations where /2N:1U — 2/<r = §96i 2<?6 and 3ga, 
respectively. Three most regular regions are indicated by A, B 
and C. 

- 45.779a; 4 + 0.230y 4 + 2.314x 2 t/ 2 (7) 

About one thousand orbits have been used to calculate the 
coefficients and the precision of the fitting is guaranteed by 
the reduced x 2 = 9.570 x 10" 16 , 3.799 x 10" 16 and 3.646 x 
10 -13 for g,s and f a respectively. 

The value 30.219 in x is assumed to be the semimajor 
axis at the center of tadpole orbits around the L5 point. 
The analytical formulas are symmetric with respect to this 
center. As we have mentioned at the end of Section 3.1, 
this symmetry is not "exactly" conserved actually. In ad- 
dition, one should note that the analytical expressions are 
more accurate and reliable in the central part of the (oo,io) 
plane where the proper frequencies are more accurately de- 
termined. Near the border between the stable and unstable 
region, the motion is more chaotic and the precision of the 
frequencies is relatively poor. 

4.3 Secular resonances map 

Given the analytical expressions of the proper frequencies, 
we can identify the main secular resonances on the initial 
plane (ao, io). For example, a Kozai resonance happens when 
g = s (since u> = w — Q, g = s infers tl) = 0), so that the posi- 
tion where the Kozai resonance may arise can be calculated 
by solving the equation g(x,y) — s(x,y) = 0. In Fig. 11, we 
plot the major secular resonances on the initial plane. The 
dynamical map is plotted as the background simultaneously, 
providing a convenient comparison between the secular res- 
onances and the dynamical stability. 

The very good agreement, between the fine structures 
in dynamical map and the locations of different secular res- 
onances, is obvious. As shown in Fig. 11, it is clear that 
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the unstable gap around io ~ 44° is determined by the v& 
resonance, the Kozai resonance can be found at high in- 
clination, and the nodal resonance i/ig takes place at low 
inclination. On one hand, around the exact location, every 
resonance has a "width" in which the resonance happens 
and the corresponding resonant argument oscillates with a 
definite amplitude. On the other hand, Fig. 11 gives only 
the location in the phase space where the secular resonance 
possibly occurs. Whether a resonance actually happens de- 
pends on whether the motion is strongly influenced by other 
mechanisms than the resonance itself. For example, the V- 
shape curves in Fig. 11 indicate the position where the vie, 
secular resonance may happen, but in fact, we only observed 
this resonance at low inclination with io < 1.5°. At higher 
inclination this resonance is hidden by other mechanisms. 

From the dynamical spectrum of g in Fig. 9, we have 
found the equality between the frequencies of f2N-.1v — 2/ CT 
and ga plays a role in forming the arc structure in the dy- 
namical map. The curve f2N:iu — 2/<r = ga in Fig. 11 co- 
incides very well with the stable arc, while curves f2N-.1v — 
2/ CT = go — g and f2N-.1v — 2j " a = ga+g, which are not plotted 
explicitly, coincide very well with the unstable arcs by the 
inner and outer sides of the stable one. Further calculations 
show that different commensurabilities between f2N-.1v — 2/ CT 
and go are most probably responsible for the multiplet arc 
structures, as other three curves labeled 1, 2 and 3 show 
in Fig. 11. Even the less regular "holes" at low inclination 
(io < 5°) around ao = 30.13 and 30.31 AU seem to be re- 
lated with this resonance family. Our calculations show that 
f2N-.1v — 2/ CT = l(/6 and s = gs assemble in these neigh- 
bourhood. One may note that these commensurabilities are 
not resonances since the d'Alembert rule of invariance to 
rotations is not fulfilled. This is because we cannot identify 
the contributions of all the longitude precessions, especially 
those with lowest frequencies. However, the plot of these 
commensurabilities gives good estimate of the position of 
the real resonances. 

As shown above, Uranus may put its influence on the 
stability of Neptune Trojans through the quasi 2:1 mean 
motion resonance. Similar effects of the quasi 5:2 resonance 
between Jupiter and Saturn (the Great Inequality) on the 
ma in-belt asteroids in resonance have been dis cussed e.g., 
in (|Ferraz-Mello. Nesvornv fc Michtchen ko 1998). We argue 
that Saturn is also very important because its apsidal pre- 
cession go is deeply involved in shaping the stable region. 

Surely we could not plot all the possible resonances 
(commensurabilities). However, we plot 2s = sg in Fig. 10, 
because it defines the upper limit io ~ 35° of stable region 
B. In a typical power spectrum of cos a, k or q, we can detect 
some less important peaks and they may evolve into some 
less important secular resonances. We list a few of these pos- 
sible commensurabilities: g = gs — 37, 2g = g 5 , 2g — gr, 2g — 
9&-,9 = 97 — 2<?8, 2s = S7 — 3s8, 3s = S7 — 2ss and 4s = S7 — s$. 
Among them, many cross the area with io G (12°, 22°) be- 
tween the stable region A and B. They themselves or their 
overlapping between each other bring irregularity to the mo- 
tion and therefore this area looks less stable than the neigh- 
boring region A and B. 



5 CONCLUSION 

Since the first Neptune Trojan was found in 2001 their num- 
ber steadily increases and now we have knowledge of 6 such 
asteroids around the Lagrange point L4. It is interesting to 
note that two of them are on highly inclined orbits. Hence we 
study in this paper the orbital stability of Neptune Trojans, 
with special interests on the inclined orbits. 

Wc first verified the symmetry between the L4 and L5 
points. We found that orbits around these two points have 
the same stability. The only difference between them is in 
the value of the osculating semimajor axis of orbits at the 
libration center in the Trojan clouds around the L4 and L5 
points. This difference was found due to the asymmetrical 
selection of initial conditions. To clarify this symmetry is 
important, not only because some papers argued that the L4 
and L5 are dynamically asymmetrical to each other, but also 
because a specific formation history of Trojan clouds may 
affect the symmetry property. If future observing confirms 
the symmetry or asymmetry, it will put strong constrains on 
the formation scenario, which is tightly related to the early 
dynamical evolution of the outer solar system. 

Using the spectral number as an indicator, we portrayed 
the dynamical map on the initial plane (ao,io)- We found 
that the inclination of stable orbits could be as high as 60° . 
Three most stable areas were located in the dynamical map, 
region A with i G (0°,12°), B with i G (22°, 36°) and C 
with io G (51°, 59°), where the amplitude of the oscillating 
resonant argument is in the range of (20° — 60°). In these 
regions more Trojan objects may be observed in future. The 
region A and B are connected to each other by an area of less 
regular orbits, but region C is isolated from A and B by an 
unstable gap at io ~ 44° . We argue that any Trojan found in 
region C must be a relic of the primordial Trojan cloud. At 
the median inclination 12° — 22°, orbits are less stable but 
many of them can survive in the outer solar system up to 
10 9 yr, and this area may host less Trojans than the stable 
regions. Trojans in this area would wander around for very 
long time before they leave the resonance. 

By calculating the dynamical spectra and the proper 
frequencies of the resonance argument, the apsidal preces- 
sion and the nodal precession, we figured out the secular 
resonances that generate the fine structures in the dynami- 
cal map. The upper limit of stable space in inclination about 
io = 60° was found to be caused by the Kozai resonance, 
in which the perihelion argument librates around uj ~ 320° , 
other than the well-known values of 0°(180°) or 90° (270°). 
The vg secular resonance was identified in the unstable gap 
around io ~ 44°, where the difference between the longi- 
tudes of the Trojan and Neptune zu — w& oscillates around 
~ 270°. We also found that the nodal secular resonance 
vis only takes place at low inclination. By constructing the 
semi-analytical expressions of the proper frequencies, we an- 
alyzed the delicate structures seen in the dynamical map on 
the (ao, io) plane. Particularly, we found that the commen- 
surabilities between the frequencies of the quasi 2:1 mean 
motion resonance /2Ni:U, of the 1:1 resonant argument f a 
and of the apsidal precession of Saturn go, are responsible 
for the multiplet arc structures in the dynamical map. 

We also check how the stability of a Trojan's orbit de- 
pends on the initial eccentricity. The preliminary results 
show that for most inclination values, the orbit needs a small 
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initial eccentricity to be stable. The largest initial eccentric- 
ity that a stable orbit may endure is eo ~ 0.10, 0.12 and 0.03 
at initial inclination around 10°, 35° and 55° respectively in 
the three stable areas (region A, B and C). A through anal- 
ysis on the eccentric orbits will be presented in the coming 
paper. 
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